X <- rnorm(1000, 1, 0.5)
Y <- rnorm(1000, 10, 2)
Z <- rbind(X,Y)
hist(Z, breaks=30, xlab="Value")
set.seed(seed=2)
X1 <- rnorm(10000, 1, 0.3)
X2 <- rnorm(10000, 10, 0.6)
trans_probs <- matrix(data=c(0.95, 0.05, 0.01, 0.99),byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;
## Load data into DataFrame
data <- data.frame(X1, X2)
n = 10000
t_states <- s0;
for (i in 1:(n-1)){
s <- sample(states, size=1, prob=trans_probs[s, ])
t_states <- c(t_states, s)
}
traj <- sapply(1:n, function(x) data[x, t_states[x]])
plot(traj[1:1000], pch=16, ylab="Value", xlab="Time")
plot(1:1000, traj[1:1000], type='l', main='Sample trajectory', ylab=c("Value"), xlab=c("Time"))
plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', ylab=c("Value"), xlab=c("Time"))
plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", ylab=c("Value"), xlab=c("Time"))
X <- rnorm(10000, 1, 0.5)
Y <- rnorm(10000, 2, 0.6)
Z <- rbind(X,Y)
hist(Z, breaks=30, xlab="Value")
##Establish Two-State Time Course Data
set.seed(seed=3)
X1 <- rnorm(10000, 1, 0.5)
X2 <- rnorm(10000, 2, 0.5)
trans_probs <- matrix(data=c(0.95, 0.05, 0.01, 0.99),byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;
## Load data into DataFrame
data <- data.frame(X1, X2)
n = 10000
t_states <- s0;
for (i in 1:(n-1)){
s <- sample(states, size=1, prob=trans_probs[s, ])
t_states <- c(t_states, s)
}
traj <- sapply(1:n, function(x) data[x, t_states[x]])
plot(traj[1:1000], pch=16, xlab="Time", ylab="Value")
plot(1:1000, traj[1:1000], type='l', main='Sample trajectory', xlab="Time", ylab="Value")
plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Time", ylab="State")
plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Time", ylab="State")
##Use 1D measures so that they are gaussian
x <- Sox2_data_unsync_20sec
L_fitX <- fitdistr(x[x$Cell_Line == "121316.2E1_ESC", "Corrected X_Distance (um)"], densfun = "normal")
L_fitY <- fitdistr(x[x$Cell_Line == "121316.2E1_ESC", "Corrected Y_Distance (um)"], densfun = "normal")
step <- seq(-1,1,0.01)
hist(x[x$Cell_Line == "121316.2E1_ESC", "Corrected X_Distance (um)"], freq = FALSE, ylim=c(0,5), xlab="X Distance (um)", xlim=c(-0.5,0.5), main="Looped Model X Distance Distribution")
L_model <- dnorm(step, mean=L_fitX$estimate["mean"], sd=L_fitX$estimate["sd"])
lines(step, L_model, lty="dashed", col="red", lwd=2)
hist(x[x$Cell_Line == "121316.2E1_ESC", "Corrected Y_Distance (um)"], freq = FALSE, ylim=c(0,5), xlab="Y Distance (um)", xlim=c(-0.5,0.5), main="Looped Model Y Distance Distribution")
L_model <- dnorm(step, mean=L_fitY$estimate["mean"], sd=L_fitY$estimate["sd"])
lines(step, L_model, lty="dashed", col="red", lwd=2)
U_fitX <- fitdistr(x[x$Cell_Line == "013116.2G8_ESC", "Corrected X_Distance (um)"], densfun = "normal")
U_fitY <- fitdistr(x[x$Cell_Line == "013116.2G8_ESC", "Corrected Y_Distance (um)"], densfun = "normal")
step <- seq(-1,1,0.01)
hist(x[x$Cell_Line == "013116.2G8_ESC", "Corrected X_Distance (um)"], freq = FALSE, ylim=c(0,2.5), xlab="X Distance (um)", xlim=c(-0.75,0.75), main="Unlooped Model X Distance Distribution")
U_model <- dnorm(step, mean=U_fitX$estimate["mean"], sd=U_fitX$estimate["sd"])
lines(step, U_model, lty="dashed", col="red", lwd=2)
legend(x="topright", lwd=2, lty="dashed", legend="Fit Gaussian", col="red")
hist(x[x$Cell_Line == "013116.2G8_ESC", "Corrected Y_Distance (um)"], freq = FALSE, ylim=c(0,2.5), xlab="Y Distance (um)", xlim=c(-0.75,0.75), main="Unlooped Model Y Distance Distribution")
U_model <- dnorm(step, mean=U_fitY$estimate["mean"], sd=U_fitY$estimate["sd"])
lines(step, U_model, lty="dashed", col="red", lwd=2)
legend(x="topright", lwd=2, lty="dashed", legend="Fit Gaussian", col="red")
##Establish as hypothetical two-state time course based on the data
set.seed(seed=4)
L_X <- rnorm(10000, mean=L_fitX$estimate["mean"], sd=L_fitX$estimate["sd"])
U_X <- rnorm(10000, mean=U_fitX$estimate["mean"], sd=U_fitX$estimate["sd"])
L_Y <- rnorm(10000, mean=L_fitY$estimate["mean"], sd=L_fitY$estimate["sd"])
U_Y <- rnorm(10000, mean=U_fitY$estimate["mean"], sd=U_fitY$estimate["sd"])
trans_probs <- matrix(data=c(0.95, 0.01, 0.01, 0.99),byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;
## Load data into DataFrame
Xdata <- data.frame(L_X, U_X)
Ydata <- data.frame(L_Y, U_Y)
n = 10000
t_states <- s0;
for (i in 1:(n-1)){
s <- sample(states, size=1, prob=trans_probs[s, ])
t_states <- c(t_states, s)
}
traj_X <- sapply(1:n, function(x) Xdata[x, t_states[x]])
traj_Y <- sapply(1:n, function(x) Ydata[x, t_states[x]])
plot(traj_X[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="X Distance (um)")
plot(traj_Y[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="Y Distance (um)")
plot(1:1000, sqrt((traj_X[1:1000])^2 + (traj_Y[1:1000])^2), type='l', main='Sample trajectory', xlab="Frames", ylab="XY Distance (um)")
plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Frames", ylab="State")
plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Frames", ylab="State")
trans_probs <- matrix(data=c(0.99, 0.01, 0.01, 0.99),byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;
## Load data into DataFrame
data <- data.frame(X, Y)
n = 10000
t_states <- s0;
for (i in 1:(n-1)){
s <- sample(states, size=1, prob=trans_probs[s, ])
t_states <- c(t_states, s)
}
traj <- sapply(1:n, function(x) data[x, t_states[x]])
plot(traj[1:1000], pch=16, xlab="Time", ylab="Value", main="Autocorrelated time course dataset derived from two independent states")
plot(1:1000, traj[1:1000], type='l', main='Sample trajectory', xlab="Time", ylab="Value")
plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Time", ylab="State")
plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Time", ylab="State")
x <- Sox2_data_unsync_20sec
images <- unique(x$C1_Image)
for (i in 1:length(images)) {
start = 1;
end = 100;
cur_image = images[i];
cur_data = x[x$C1_Image == cur_image,];
for (j in start:end) {
if (nrow(cur_data[cur_data$`C1_T Value (frame)` == (j - 1),]) != 0 & nrow(cur_data[cur_data$`C1_T Value (frame)` == j,]) != 0){
x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_X.Displacement..um."] <- (x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == (j - 1), "Relative_C1_Corrected_X Value (um)"] - x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_C1_Corrected_X Value (um)"]);
x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_Y.Displacement..um."] <- (x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == (j - 1), "Relative_C1_Corrected_Y Value (um)"] - x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_C1_Corrected_Y Value (um)"]);
}
else{
x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_X.Displacement..um."] <- NA;
x[x$C1_Image == cur_image & x$`C1_T Value (frame)` == j, "Relative_Y.Displacement..um."] <- NA;
}
}
}
hist(x[x$Cell_Line == "101515.C3_ESC", "Relative_X.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))
hist(x[x$Cell_Line == "101515.C3_ESC", "Relative_Y.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))
hist(x[x$Cell_Line == "013116.2G8_ESC", "Relative_X.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))
hist(x[x$Cell_Line == "013116.2G8_ESC", "Relative_Y.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))
hist(x[x$Cell_Line == "121316.2E1_ESC", "Relative_X.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))
hist(x[x$Cell_Line == "121316.2E1_ESC", "Relative_Y.Displacement..um."], breaks=c(-10, seq(-0.95, 0.95, by=0.05), 10), xlim=c(-1, 1))
##Use 1D measures so that they are gaussian
L_fitX <- fitdistr(na.omit(x[x$Cell_Line == "121316.2E1_ESC", "Relative_X.Displacement..um."]), densfun = "normal" )
L_fitY <- fitdistr(na.omit(x[x$Cell_Line == "121316.2E1_ESC", "Relative_Y.Displacement..um."]), densfun = "normal")
step <- seq(-1,1,0.01)
hist(x[x$Cell_Line == "121316.2E1_ESC", "Relative_X.Displacement..um."], freq = FALSE, ylim=c(0,5), xlab="X Displacement (um)", xlim=c(-0.5,0.5), main="Looped Model X Displacement Distribution")
L_model <- dnorm(step, mean=L_fitX$estimate["mean"], sd=L_fitX$estimate["sd"])
lines(step, L_model, lty="dashed", col="red", lwd=2)
hist(x[x$Cell_Line == "121316.2E1_ESC", "Relative_Y.Displacement..um."], freq = FALSE, ylim=c(0,5), xlab="Y Displacement (um)", xlim=c(-0.5,0.5), main="Looped Model Y Displacement Distribution")
L_model <- dnorm(step, mean=L_fitY$estimate["mean"], sd=L_fitY$estimate["sd"])
lines(step, L_model, lty="dashed", col="red", lwd=2)
U_fitX <- fitdistr(na.omit(x[x$Cell_Line == "013116.2G8_ESC", "Relative_X.Displacement..um."]), densfun = "normal")
U_fitY <- fitdistr(na.omit(x[x$Cell_Line == "013116.2G8_ESC", "Relative_Y.Displacement..um."]), densfun = "normal")
step <- seq(-1,1,0.01)
hist(x[x$Cell_Line == "013116.2G8_ESC", "Relative_X.Displacement..um."], freq = FALSE, ylim=c(0,2.5), xlab="X Displacement (um)", xlim=c(-0.75,0.75), main="Unlooped Model X Displacement Distribution")
U_model <- dnorm(step, mean=U_fitX$estimate["mean"], sd=U_fitX$estimate["sd"])
lines(step, U_model, lty="dashed", col="red", lwd=2)
legend(x="topright", lwd=2, lty="dashed", legend="Fit Gaussian", col="red")
hist(x[x$Cell_Line == "013116.2G8_ESC", "Relative_Y.Displacement..um."], freq = FALSE, ylim=c(0,2.5), xlab="Y Displacement (um)", xlim=c(-0.75,0.75), main="Unlooped Model Y Displacement Distribution")
U_model <- dnorm(step, mean=U_fitY$estimate["mean"], sd=U_fitY$estimate["sd"])
lines(step, U_model, lty="dashed", col="red", lwd=2)
legend(x="topright", lwd=2, lty="dashed", legend="Fit Gaussian", col="red")
###Establish as hypothetical two-state time course based on the data
#set.seed(seed=4)
L_X <- rnorm(10000, mean=L_fitX$estimate["mean"], sd=L_fitX$estimate["sd"])
U_X <- rnorm(10000, mean=U_fitX$estimate["mean"], sd=U_fitX$estimate["sd"])
L_Y <- rnorm(10000, mean=L_fitY$estimate["mean"], sd=L_fitY$estimate["sd"])
U_Y <- rnorm(10000, mean=U_fitY$estimate["mean"], sd=U_fitY$estimate["sd"])
trans_probs <- matrix(data=c(0.98, 0.02, 0.02, 0.98), byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;
## Load data into DataFrame
Xdata <- data.frame(L_X, U_X)
Ydata <- data.frame(L_Y, U_Y)
n = 10000
t_states <- s0;
for (i in 1:(n-1)){
s <- sample(states, size=1, prob=trans_probs[s, ])
t_states <- c(t_states, s)
}
traj_X <- sapply(1:n, function(x) Xdata[x, t_states[x]])
traj_Y <- sapply(1:n, function(x) Ydata[x, t_states[x]])
plot(traj_X[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="X Displacement (um)")
plot(traj_Y[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="Y Displacement (um)")
plot(1:1000, sqrt((traj_X[1:1000])^2 + (traj_Y[1:1000])^2), type='l', main='Sample trajectory', xlab="Frames", ylab="XY Distance (um)")
plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Frames", ylab="State")
plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Frames", ylab="State")
#set.seed(seed=4)
L_X <- x[x$Cell_Line == "121316.2E1_ESC", "Relative_X.Displacement..um."]
U_X <- x[x$Cell_Line == "013116.2G8_ESC", "Relative_X.Displacement..um."]
L_Y <- x[x$Cell_Line == "121316.2E1_ESC", "Relative_Y.Displacement..um."]
U_Y <- x[x$Cell_Line == "013116.2G8_ESC", "Relative_Y.Displacement..um."]
trans_probs <- matrix(data=c(0.98, 0.02, 0.02, 0.98), byrow=TRUE, nrow=2)
states <- c(1, 2)
s0 <- 1;
s <- s0;
## Load data into DataFrame
Xdata <- data.frame(L_X[1:5000], U_X[1:5000])
Ydata <- data.frame(L_Y[1:5000], U_Y[1:5000])
n = 5000
t_states <- s0;
for (i in 1:(n-1)){
s <- sample(states, size=1, prob=trans_probs[s, ])
t_states <- c(t_states, s)
}
traj_X <- sapply(1:n, function(x) Xdata[x, t_states[x]])
traj_Y <- sapply(1:n, function(x) Ydata[x, t_states[x]])
plot(traj_X[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="X Displacement (um)")
plot(traj_Y[1:1000], pch=16, ylim=c(-0.8,0.8), xlab="Frames", ylab="Y Displacement (um)")
hmm <- depmix(response = list(traj_X ~ 1, traj_Y ~ 1), ntimes=5000, nstates=2, family=list(gaussian(), gaussian()), instart=c(1,0))
f <- fit(hmm)
## iteration 0 logLik: 6696.141
## iteration 5 logLik: 7041.577
## iteration 10 logLik: 7166.998
## iteration 15 logLik: 7228.651
## iteration 20 logLik: 7250.898
## iteration 25 logLik: 7257.632
## iteration 30 logLik: 7259.593
## iteration 35 logLik: 7260.171
## iteration 40 logLik: 7260.346
## iteration 45 logLik: 7260.4
## iteration 50 logLik: 7260.418
## iteration 55 logLik: 7260.423
## iteration 60 logLik: 7260.425
## converged at iteration 65 with logLik: 7260.425
summary(f)
## Initial state probabilties model
## pr1 pr2
## 1 0
##
## Transition matrix
## toS1 toS2
## fromS1 0.954 0.046
## fromS2 0.121 0.879
##
## Response parameters
## Resp 1 : gaussian
## Resp 2 : gaussian
## Re1.(Intercept) Re1.sd Re2.(Intercept) Re2.sd
## St1 0 0.089 -0.001 0.090
## St2 0 0.175 0.000 0.178
esttrans <- posterior(f)
plot(1:1000, sqrt((traj_X[1:1000])^2 + (traj_Y[1:1000])^2), type='l', main='Sample trajectory', xlab="Frames", ylab="XY Distance (um)")
plot(1:1000, t_states[1:1000], type='l', col="red", main='Actual State', xlab="Frames", ylab="State")
plot(1:1000, esttrans[1:1000,1], type="l", col="blue", main="Predicted State", xlab="Frames", ylab="State")